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Abstract An error analysis result is given for classical Gram-Schmidt fac- 
torization of a full rank matrix A into A = QR where Q is left orthogonal 
(has orthonormal columns) and R is upper triangular. The work presented 
here shows that the computed R satisfies R T R = A T A + E where E is an 
appropriately small backward error, but only if the diagonals of R are com- 
puted in a manner similar to Cholesky factorization of the normal equations 
CNJ ■ matrix. 

A similar result is stated in [Giraud at al, Numer. Math. 101(1):87- 
100,2005]. However, for that result to hold, the diagonals of R must be com- 
puted in the manner recommended in this work. 



OO 

in 

<N 
\o . 

The classical Gram-Schmidt (CGS) orthogonal factorization is analyzed 
■ in a recent work of Giraud et al. [5] and in a number of other sources [UHJ 

O ; EEEE], US §6.9], §2.4.5]. 

r"| ■ For a matrix A £ ]R mxn ( m > n ) with rank(A) = n, in exact arithmetic, 

' the algorithm produces a factorization 

£■ A = QR (1) 

> 

where Q is left orthogonal (i.e. Q Q = I n ), and R £ M. nxn is upper trian- 
' gular and nonsingular. In describing the algorithms, we use the notational 
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conventions, 

A = (a 1: . . . ,a n ), Q = (qi, . . . ,q„), 
R = (r jk ). 

The algorithm forms Q and R from A column by column as described 
in the following pseudo-code. We label this algorithm CGS-S, for classical 
Gram-Schmidt "standard." 

Algorithm 1 (Classical Gram Schmidt Orthogonal Factorization (Standard) (CGS-S)) 

rn = ||ai|| 2 ;qi = ai/ru; 
Ri = (rn);Qi = (qi); 
for k = 2 : n 

Vfe = afc — Qfe-iSfe; 
rkk = ||vfc|| 2 ; 
qfe = v fe /r fefe ; 

fc- 1 1 

p _ k-l f R k -i s fc \ 

end; 

Q Qni R Rn: 

As is well known [2j p. 63, §2. 4. 5], in floating point arithmetic, Q is far from 
left orthogonal. The authors of [5] prove a number of results about classical 
Gram-Schmidt. This note shows that for one of their results (Lemma 1 in [5]), 
the diagonal elements r^k should be computed differently from Algorithm [TJ 
substituting a Cholesky-like formula for rkk rather than setting rkk = \\^k lb- 
That change produces the Algorithm [2] Since it uses a Pythagorean identity 
to compute the diagonals of R, we call it CGS-P for "classical Gram-Schmidt 
Pythagorean." . 

Algorithm 2 (Cholesky like Classical Gram Schmidt Orthogonal Factorization (CGS-P)) 

ni = || a i|| 2 ;qi = ai/m; 
Ri = (ru);Qi = (qi); 
for k = 2 : n 

s*; = Ql-i&k; 
Vfe = afe - Qk-iSk] 
■0fc = 1 1 a fc j 1 2 ; 4>k = 1 1 s fc j 1 2 ,- 

Tkk = (^-0fc) 1/2 (ipk+lpk) 1 ' 2 

Qk = v fe /r fefe ; 

k-l 1 

k-l f R k -i s fe \ 

^ = 1 ^ r fefe J> 0fe 



fc- 1 1 
= (Qfe-i qfe); 



k-l 1 

= (Qk-i qfe); 
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end; 

Q — Qn 5 R Rn i 

We assume that we are using a floating point arithmetic that satisfies the 
IEEE floating point standard. In IEEE arithmetic 

f£(x + y) = (x + y)(l+6), \6\<e M 

for results in the normalized range [H p. 32]. 

Letting em be the machine unit, we follow Golub and Van Loan [6l §2.4.6] 
and use the linear approximation 

(l + e M ) p{n) = l+p{n)e M + 0{e 2 M ) 

for a modest function p(n) thereby assuming that the 0(e 2 m ) makes no sig- 
nificant contribution. 

For the sake of self containment, we give Lemma 1 from [5]. 

Lemma 1 In floating point arithmetic with machine unit Em, the com- 
puted upper triangular factor from Algorithm]]^ satisfies 

R T R = A T A + E, \\E\\ 2 < c(m,n)\\A\\le M 

where c(m,n) = 0(mn 2 ). 

As stated, this lemma is not correct for Algorithm [TJ but a slightly dif- 
ferent version of this result holds for Algorithm [2j 
We define the four functions 

1 k = 1 

2V2mk + 2y/k k = 2, . . . , n, 

( , \ / m + 2 k = 1 , . 

C2[m > fcj - \ 3.5mk 2 - 1.5mfc + 16fc fc = 2, . . . [Z) 

c^(m, k) = 0.5c 2 (m, fc), 04(771, fc) = c 2 (m, fc) + 2ci(m, fc), 
we let Ak be the first fc columns of A, and let 

K 2 {R k ) = \\Rkh\\Rk%- 

The new version of Lemma [T] is Theorem [1] 

Theorem 1 Assume that in floating point arithmetic with machine unit em , 
for the R resulting from Algorithm^ for each k, we have 

c 4 {m,k)£ M K 2 {R k ) 2 < 1. (3) 

Let Ak G R mxfc consist of the first fc columns of A. Then, for fc = 1, . . . , n, 
to within terms of 0{e 2 m ), the computed matrices Rk and Qk satisfy 

QkRk — Ak = AA kl \\AAkh < ci (m, A:)||Afc|| 2 e A/ , (4) 

Rl R k - A\ A k = E k , \\Ek\\ 2 < c 2 {m,k)\\Akh 2 eM, (5) 

\\Rkh = \\AhhQ- +Mfc). iMfcl < c 3 (m,k)£ M , (6) 

\\I-QlQkh < c 4 (m,k)K 2 {R k ) 2 EM, (7) 

IIQfclla < V2. (8) 



Ci(m, fc) 
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The proof of Theorem [T] is given in the appendix. 

The restriction (|3|) assures that R is nonsingular, and that (|7|) and ([8]) 
hold. A weaker assumption that assures that R is nonsingular and that ||Q/c||2 
is bounded would yield bounds similar to (|4]), ([5]), and (J6|). 

Remark 1 The condition ([3|) and the bound ([7|) are stated in terms of K 2 (R k ). 
We now show how it may be stated in terms of 

K 2 (A k ) = \\A k \\ 2 \\Al\\ 2 

where A\, is the Moore-Penrose pseudoinverse of A k . In exact arithmetic, 
K 2 {A k ) and n 2 (R k ) are the same quantity, and equation (J6j) states that ||-Rfc||2 
and ||Afc||2 are nearly interchangablc in floating point arithmetic. To relate 
\\R~k 1 ||2 and ll^fclU, we use eigenvalue inequalities. 
From the fact that 

HflkV = V^C^)' \\A\\\^ = ^\ k {AlA k ) (9) 

where X k (-) denotes kth largest (and therefore smallest) eigenvalue, we can 
obtain an upper bound for ||^4j.||2 using Weyl's monotonicity theorem [101 
Theorem 10.3.1]. Applying that theorem to ([5]), we have 

Xk(RlRk) > X k (Al A k ) - \\E k \\ 2 

> \ k {Al A k ) - e M c 2 {m, k)\\A k g + 0(e 2 M ) 
= \ k (Al A k ) - e M c 2 (m, k)\\R k f 2 + 0{e 2 M ) 

> X k (Al A k )(l - ( k ) 

where 

Ck - e M c 2 (m, k)n 2 {R k f + 0{e 2 M ). (10) 

Using ([9]), we have 

\\Rlh<\\A- k %{l-Ck)- 1/2 . 
From (|6|), we may conclude that 

K 2 (R k ) < K 2 (A k )(l + fi k )(l - ( k y 1/2 . 

Thus a slight variation of the condition (J3j> may be stated in terms of K 2 {A k ). 
Since it fits more naturally into the proof of Theorem [T] and it is more easily 
computed than n 2 (A k ), we use n 2 (R k ). 

The conclusion of Theorem [1] does not hold for Algorithm [TJ as shown by 
the following example. We were able to construct several similar examples. 
Both examples were done in MATLAB version 7 on a Dell Precision 370 
workstation running Linux. 

Example 1 We produced a 6 x 5 matrix with the following MATLAB code. 
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Algorithm 


\\A*A-lCR\\ a /\\A\\i 


\\I-Q' 1 Qh 


CGS-S (Algorithm Hi 


4.5460e-9 


3.9874e-6 


CGS-P (Algorithm |2| 


3.3760e-17 


5.2234e-5 



Table 1 Orthogonality and Normal Equations Error from CGS Algorithms for 
Example [T] 



B=hilb(6); 

Al = ones(6, 3) + B{ : , 1 : 3) * le - 2; 
B=pascal(6); 
A2 = B(:,1: 2); 
A=[A1 A2]; 

The command hilb(6) produces the 6x6 Hilbert matrix, the command 
ones(6,3) produces a 6 x 3 matrix of ones, and the command pascal(6) pro- 
duces a 6 x 6 matrix from Pascal's triangle. The condition number of R from 
Algorithm 2,k 2 (R) = ||i?|| 2 ||.R _1 ||2, computed by the MATLAB command 
cond, is 3.9874 • 10 6 , thus given that e M ~ 2.2206 • 1CT 16 in IEEE double 
precision, R is neither well-conditioned nor near singular. 

We computed the Q-R factorization using Algorithm [1] (CGS-S) and 
then we computed the same factorization using Algorithm [2] (CGS-P). The 
resulting Q and R satisfy the results in Table [T] 

The bound on \\A T A — R T R\\% in appears to be satisfied if rkk is 
computed as in Algorithm [2 but it is not if rkk is computed as in Algorithm 

m 

A larger, more complex, but better conditioned example is given next. 

Example 2 A large class of examples where CGS-S obtains a large value of 
\\A T A - R T R\\ 2 /{\\A\\l), but CGS-P arises from glued matrices. A general 
MATLAB code for these glued matrices is given by 

function [A] =create_gluedmatrix (condA_glob , condA,m,nglued,nbglued) 
n = nglued*nbglued; 
A = orth(rand(m,n) ) ; 

A = A*diag([10.~(0:condA_glob/(n-l) : condA_glob)] ) *orth(randn(n,n) ) ; 
ibeg = 1; 
iend = nglued; 
for i=l:nbglued, 

A( : , ibeg: iend) = A( : , ibeg : iend) *diag( [10 . ~ (0 : condA/ (nglued- 1) : condA)] ) . . . 

*orth (randn (nglued, nglued) ) ; 
ibeg = ibeg+nglued; 
iend = iend+nglued; 
end 

Here m represents the number of rows of A, nglued is the number of columns 
in a block, nbglued is the number of blocks that are glued together, and 
n = nglued x nbglued is the number of columns in the matrix. The parameter 
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Algorithm 


WA 1 ' A-R 1 ' R\\ 2 /\\A\\i 


\\i-Q r Qh 


CGS-S (Algorithm Hi 


3.8744e-6 


9.3676e-4 


CGS-P (Algorithmic 


2.8729e-16 


1.8972e-12 



Table 2 Orthogonality and Normal Equations Error from CGS Algorithms for 
Example [2] 



condA is the condition number of a block, and condA_glob is a parameter to 
couple the blocks together. The MATLAB command orth(X) produces an 
orthonormal basis for the range of X , thus the command ortri(raridn(m,n) ) 
produces a random orthogonal matrix. 

For this example, we used the parameters 

condA_glob = 1; condA = 2; m = 200; nglued = 5; nbglued = 40; 

for which we obtained a 200 x 200 matrix with condition number 506.92 (the 
condition number of the orthogonal factor R is about the same). We also used 
the command randn( ' state ' ,0) to reset the random number generator to 
its initial state. Table [5] summarizes the results from applying CGS-S and 
CGS-P to this matrix. 

For this example, the loss of orthogonality of CGS-S is far in excess of 
0(en 2 (R) 2 ), whereas the loss of orthogonality for CGS-P is well within that 
bound. The error \\A T A — R T R\\ 2 is far larger for CGS-S than it is for 
CGS-P and is much greater than 0(£m||-A||2)- 



Conclusion 

The upper triangular factor R from classical Gram-Schmidt has been shown 
to satisfy the bound ([5|) provided that the diagonal elements of R are com- 
puted as they are in the Cholesky factorization of the normal equations ma- 
trix. If these diagonal elements are computed as in standard versions of clas- 
sical Gram-Schmidt, no bounds such as ([5]) or ([7]) may be guaranteed. 
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Appendix. Proof of Theorem Q] 

To set up the proof of Theorem [TJ we require a lemma. 

Lemma 1 Let Q € R mxn and R e R nxn be the results of Algorithm^ in 
floating point arithmetic with machine unit Em and that R satisfies ^ . Then 

r u = ||ai|| 2 (l + <5i), \5i\<(0.5m+l)s u + O(e 2 M ) (11) 
and for k = 2, . . . , n 

r kk = (||a fe ||2(l + 4)-||s fe ||2(l + zA fe )) 1/2 , (12) 
\5 k l\A k \ < ( m + 8)e M + 0(e 2 M ), 

\\s k h< ||a fe || 2 (l + C), ICI < (m + 2)e M + 0(e 2 M ). (13) 

Proof. Equation (fTTj) is just the error in the computation of | j ax 1 1 2 - In the 
computation of r kk , k = 2, . . . , n, note that 

^=MI|afc|| 2 ) = ||a fc || 2 (l + e[* ) ) > (14) 

& = /mM 2 ) = IM2(i + 4 fc) ), (is) 

|ej fc) | < (0.5m+l)e M + O(e 2 M ), » = 1,2. 

Using ([3]), we conclude that R is nonsingular, thus r kk > for all fc. Thus 
in Algorithm^ r kk > only if vb k > 4>k- 
To get ([12]). note that 

r kk = v%-<W# + <Ml + 4 fe) ), l4 fe) | < 3£M + 0(e 2 M ). 
Thus using (fTl)) and (fll))) . we have 

rkk = V !|a fc !||(l + £l fc) ) 2 - !|s fc !||(l + e ( fc) )2(l + 4 fe )) 
= (\\ ak \\ 2 (l + 5 k )-\\s k \\ 2 (l + A k )) 1/2 

where 

4 = (i + 4 fc) ) 2 (i + 4 fc) ) 2 -i ; 

4 fe = (l + e f) 2 (l + ef) 2 -l. 
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That yields 

\S k \,\A k \ < {m + S)e M + 0(e 2 M ). 

Therefore r kk satisfies (fT2|). 

Since t/j k > 4>k as outlined above, from (fT4|) -(ri5 |) . we have 

A = ||a fe || 2 (l + e[ k) ) > cp k = ||s fc || 2 (l + 4 fe) ) 

thus 

n Sfc n 2 < na fc || 2 (i + e t fc >)(i + 4 fc) )- 1 

< ||afc|| 2 (l + 

where £ satisfies (|T3")) . □ 

As a consequence of the singular value version of the Cauchy interlace 
theorem [BJ p.449-450, Corollary 8.6.3], we have that ||i? fe || 2 < ||-R|| 2 and 
ll-Rfc 1 |b < || R 1 1| 2- We will use these facts freely in the proof of Theorem [TJ 

We can now prove Theorem [T] 

Proof, [of Theorem [T] The results ((I])-© are proven by induction on k. 
First, consider k = 1. From Lemma [TJ we have ([TTj) . so 

r n = ||ai||a(l + *i), \Si\ < (0.5m + l)e M + O{e 2 M ) 
which implies that 

RlRi=r 2 n = ||a 1 ||2(l+^ 1 ) 2 

= AlA 1 {l + 5 l f=AlA 1 +E 1 

where 

Ei = 2SxAj Ax + 5lA^A x . 

Thus 

ll^ilb = \Ex\ < (m + 2)||a 1 ||| £M + 0( e | f ) = (m + 2)\\A 1 \\ 2 2 e M + 0{e 2 M ). 
Also, we can conclude from standard error bounds that 
qi = (I + Gi)ai/r Uj ||Gi|| 2 < s M - 

Therefore 

Ai - Q 1 R 1 =a 1 qim = -Giai 

so that 

\\Ai - QiRi\\ 2 = \\ai - qjriiHa < ||Gi||a||ai|| a < £j\/||ai|| 2 - (16) 

Assume that (H])-® hold for k — 1, and prove them for k. We first prove 
(HI)-©, and then show that ©-© follow. 

First, we start with error bounds of the computation of the vectors s k ,v k , 
and q k to prove (g]). Note that 

Sfc = fHQLi**) = Ql-i*k - Ss k (17) 
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where 



Also, we have 



||<5s fe || 2 < mVk- l||Qfc-i||2||a fc ||2£M + 0(e 2 M ) 
< V2(*-l)m||afc|| 2 eM + 0(e 2 M ). 

v fc = fl(a k - Qk-iSk) = a fe - Qfe-is fe - 5v fe 



(18) 
(19) 



where 



||<5v fe || 2 < ||a fc || 2 e M + Vfc^lm||Q fc _i|| 2 ||s fc || 2 e M + Q(e 2 M ). 



^From (|13[) . the bound on ||sfc|| 2 in (flUI) , and the induction hypothesis on 
Qk-i, we have 



||5v fe || 2 < (y/2(k - l)m + l)||a fc || 2 £ M + 0(4/)- 
Again using the bound on ||sfc|| 2 in fTS")) . we note that 



(20) 



|v fc + <rv fe || 2 = 



< 
< 



< 



< 



ll a *v 
ll a fc 
ll a fc 
ll a fc 
l a fr 
l a fr 
ll a fc 

ll a 'v 



l - 2a^g fe _is fe + IIQfc-iSfcUl 
2-2||s fe ||2 + ||g fe _ lSfc ||2-2( ( 5s fc ) T s fc 



2-2||s fe ||2 + ||Q fc _ 1 ||2|| Sfc ||2_2( ( 5s fc ) T s fe 



i-2||s fc ||l + 2||s fc ||i-2(<Js fc ) r s fc 



£ - 2(8s k y s k 
| + 2||5s fe || 2 ||s fe || 2 
2 + 2||5s fc || 2 ||a fc || 2 + 0( £ ^) 

|2 i 



(1 + y/2{k-l)me M ) 2 + 0(e 2 M ). 



Thus 



||v fe || 2 < ||a fc || 2 (l + (3 ^2(k - l)m)e M ) + 0(e 2 M ) = \\a k \\ 2 + 0(e M ). 
We note that 



If we let 

then 

where 

That yields 



qfc = (I+ G fe )vfc/r fc fc, ||Gfc|| 2 < e M - 
AA k = Q k R k - A k 
AA k = (AAk-i 5a k ) 

Sa k = (I + G k )v k + Qfc-iSfc - a fc , 
= Gfcv fc - 5v k . 



\\Sa k \\ 2 < ||G fc || 2 ||v fc ||2 + ||*v fc || 2 < (2y/2{k - l)m + 2)e M ||a fe || 2 + 0{e 2 M ). 
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To bound || A/4fc||2, we give a recurrence for bounding ||Z\Afc||j? in terms 
of ||j4fe||ir, then use the bound ||^4fc||_F < V^HAfelU- We show that 

||AA fc || F < ci(m, k)\\A k \\ F e M + 0(e 2 M ). 

For k = l, 

\\AAxWf = ||ai|| 3 = e M \\aLi\\ 2 = £mPi||f- 
Using properties of the Frobenius norm, 

\\AA k f F < \\AA h _4l + \\8a k \\l 

< [c\{m, k - l)||A fc _x||| + (2V2(fe-l)m + 2) 2 ||a fc j| 2 ] £ 2 M + 0(e 3 M ) 

< max{c 2 (m, k - 1), (2 y/2(k - l)m + 2) 2 }(|| A k ^ \\ 2 F + \\a k \\ 2 2 )e 2 M + 0{ 
= c 2 (m,k)\\A k \\ 2 F s 2 I + 0(s 3 M ). 

A quick induction argument yields 

cx(m, k) = 2y/2(k- l)m + 2 < 1\F2km + 2. 

Thus 

||AAfc||a < ||^A fc || F < 6i(m,fc)e J if||i4 fc ||p.+0(e5 tf ) < \/fcci(m, fc)|| A fc || 2 +0(e^) 

yielding (gj) with Ci(m, fc) = 2\f2mk + 2\[k > Vkci(m, k). 
To prove ([5]), note that 

k - 1 1 

*-**-«*-*-» (»? 2) 

where using Lemma [TJ we have 

w fc = i?Li s fc ~ Al_ x & kl 

T IT 

e kk = s fc s k + r fcfe — a k a k 
= S k a k a k - A k s k s k - 

Using the bounds on S k and A k in (jT3J) , we have 

\e k k\ < |4|||a fe ||^ + |Zi fc |||s fc || 2 

<(\6 k \ + \A k \)\\a k \\ 2 + 0(e 2 M ) 

< 2(m + 8)||a fc || 2 e M + 0(£L) 

< 2(m + 8)\\A k \\ 2 e M + 0(e 2 M ). 

Since 

s fc + Ss k = Ql_ x ak, A k -i + AA k -i = Qk-iRk-i 
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we have 

w k = Rl_xS k - Aj^ak 

= R k -xQk-i a k - Rk-i 5s k ~ ^Li a fc 

= AAl_ 1 a k -Rl_ 1 5a k . (22) 

So that || Wfe 1 1 2 has the bound 

||w fc || 2 < ||2L4 fc _i|| 2 ||a fc ||2 + ||i2fc-i||2||5sfc|| 2 + 0(4/) 

< (ci(m,fc- l)||A fc _i|| 2 ||a fc || 2 + v / 2(fc-l)m||A-i||2||a fe || 2 )e A/ 

< [2V2m{k - 1) + 2Vk - 1 + V2(*-lH||i4 fc _i|| 2 ||a*|| 2 eiif + 0{e 2 M ) 
<7m{k-l)\\A k \\le M + 0{e 2 M ) (23) 

We have that 

i (V £)■■+■(,£?)'• 

< max{||£ , fc _i|| 2 , \e kk \} + ||w fc || 2 

< [max{ c 2 (m, k - 1), 2(m + 8)} + 7m(fc - 1)]|| A fc ||^ M + 0(4r) 

< [c 2 (m, k - 1) + 2(m + 8) + 7m(fc - l)}\\A k \\ 2 2 e M + 0(e 2 M ) 
<c 2 (m,k)\\A k \\ 2 2 s M + 0(e 2 M ) (24) 

where 

k 

c 2 (m, k) = £)[2(m + 8 ) + 7m U - !)] 

j'=i 

= 3.5m(fc - l)fc + 2mk + 16k. 

Thus we have the expression for c 2 (m, k) given in equation (|2"|). 

To prove ([5])-([8]). we simply apply l[4"|)-(|5|). Equation ((6]) results from 
noting that 

\\R k \\ 2 = \\R T k R k \\ 2 = \\A T k A k + E k \\ 2 

< \\A T k A k \\ 2 + \\E k \\ 2 < (l + c 2 (m,k)e M )\\A k \\ 2 2 + 0(e 2 M ). 

Thus, 

||-R*||a < {l + c 3 (m,k)e M )\\A k \\ 2 + 0(e 2 M ) 

where 

1 + c 3 (m, k)e M + 0{e 2 M ) = y/T+ c 2 (m, k), 
that is, 03(771, k) = 0.5c 2 (m, k). Reversing the roles of R k and A k yields 

\\A k h < (1 + c 3 K fc)e M )||i2ib||2 + O(e^), 

thus we have ©. 

To get (0, we note that 

Q fe = (A fc + AA k )R- 1 
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SO that 

I-Ql Qk = Rk T {Rl Rk - {A k + AA k ) T {A k + AA^R- 1 

= R k T (E k - A T k AA k - (AA k ) T A k - (AA k ) T (AA k ))R-\ 

Thus 

\\I-Ql Qkh < llVHKll^lh + mAkh\\A k \\ 2 + \\AA k \\l) 

< \\R- 1 \\l{c 2 {m 1 k)\\A k \\l + 2c 1 (to, k)\\A k \\\ + e M cl(m, k)\\A k \\ 2 2 )e M + 0(e 2 M ) 

< \\R k \\l\\R k 1 \\ 2 2(c2(rn 7 k) + 2c 1 (m,k))e M + 0(e 2 M ) 
= c A {m,k)\\R k \\l\\R~ 1 \\le M + 0{e 2 M ) 

where 04(771, k) = C2(m, k) + 2c\(m, k). 
Finally, to get ([8]), we have that 

\\Qk\\l = \\QlQkh = \\l-QlQk-lh 

< \\lh + \\l-QlQkh 
<l + \\I-QlQkh 

< 1 + c 4 (to, AOII^-ll'Pfc 'II^m + (e 2 M ) < 2 + 0(4,). 
Taking square roots yields ([5]). □ 



